Noname manuscript No. 

(will be inserted by the editor) 



Simulation and estimation for the fractional Yule 
process 

Dexter O. Cahoy • Federico Polito 



Received: 09-13-2010 / Accepted: date 



Abstract In this paper, we propose some representations of a generalized 
linear birth process called fractional Yule process (fYp). We also derive the 
probability distributions of the random birth and sojourn times. The inter- 
birth time distribution and the representations then yield algorithms on how 
to simulate sample paths of the fY p. We also attempt to estimate the model 
parameters in order for the fYp to be usable in practice. The estimation pro- 
cedure is then tested using simulated data as well. We also illustrate some 
major characteristics of fYp which will be helpful for real applications. 
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1 Introduction 

The pure birth process is undoubtedly considered as one of the simplest 
branching processes. It has a Markovian structure and has already been exten- 
sively studied in the past. When the birth rate is linear, it is then usually called 
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the pure linear birth or classical Yule or Yule-Furry process (Yp). The pure 



linear birth process has been introduced by McKendrick (1914), and has been 
widely used to model various stochastic dynamical systems such as cosmic 
showers, epidemics, and population growth to name a few. 

For the sake of completeness, we review some known properties of the 
classical Yule process which will be used in the succeeding discussion. Let %l(t) 
be the number of individuals in a Yule process with a single initial progenitor 
and birth intensity A > 0. The fcth state probability or the probability of 
having exactly k individuals p&(t) = Pr{9T(t) = k | 91(0) — 1} in a growing 
population at time t > solves the following Cauchy problem: 



f 4p fc (t) = -Afcpfe(t) + X(k - l)j> fc _i(t), k > 1, 



Pfc(0) 



1, k = l, 
0, k>l, 



(1.1) 



where po(0) = 0. The explicit solution to ( |1.1| is 

Pfc(t) = e- x \l-e- xt ) k -\ t>0, fc>l, 



with mean Em{t) = e xt . To make the Yule process more flexible in taking 
into account more complex non-Markovian behaviour, some authors (Uchaikin 



et al. (2008), Orsingher and Polito (2010)) proposed a more general model 
called the fractional Yule process (fYp). A similar generalization of other point 
processes such as the Poisson process has previously b een carried out by|Repin 
and Saichev (2000[) , | Jumarie (2001|), |Laskin (2003[), |Wang and Wen (2003f 



Mainardi et al. (2004),|Wang et al. (20061, |Wang et al. (2007), Mainardi et 



al. (2005|),|Cahoy (2007|), |Uchaikin and Sibatov (2008[ ) , |Uchaikin et al. (2008] ) 
and |Beghin and Orsingher (2009 ). 

The aim of this paper is twofold: We want to derive related representa- 
tions of fYp in terms of some classical or standard processes, and we want to 
construct algorithms on how to simulate a fYp and estimate the parameters. 

We organize the rest of the paper as follows: In Section 2, we show the 
fractional generalization of the pure linear birth process. In Section 3, it is 
illustrated that a pure linear birth process can also be viewed as a classical 
linear pure birth process with Wright-distributed random rates evaluated on 
a stretched time scale, i.e., 



W(t) = *rt s (0, ^e(o,i], 

where 5 is a random variable having the Wright probability density function 



w- v ,x- v (-z) = J2 



(-0! 

r!r(l - v(r + l))' 



(1.2) 



Furthermore, some Poisson-related representations are proved. In Section 4, 
we derive the birth and inter-birth time distributions. The structural represen- 
tation, fractional moments of the sojourn and birth times are also shown. In 
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Section 5, we generate sample paths of a fYp using our algorithms. In Section 6, 
an estimation procedure is proposed using the moments of the log-transformed 
data, and some empirical results are showed as well. Section 7 concludes the 
paper with a discussion on the key points and possible extensions of this study. 



2 Generalization of the Yule process 



The fractional generalization of the Cauchy problem (1.1 1 was first carried 
out in Uchaikin et al. (2008), Section 8, and is described as follows: The au- 



thors defined the following difference-differential equations governing the state 
probabilities p%(t) = Pr {N v (t) = k \ N"(0) = 1}: 



Pk(t) = A 



fc-i 



i=i 



r(i - v 



■5 k , u ye (0,1], k> 1, 
(2.1) 



where the initial condition 



P*(0) 



1, fc = l, 
0, k > 1, 



is incorporated into equation (2.1 ) through the Kronecker delta 5/. The frac- 



tional derivative appearing in (2.1) is the so-called Ricmann-Liouville opera- 
tor, and is defined as 



fit) 



r(l-v) dt JO (t-s) 



ds, v e (0, 1), 
u=l. 



(2.2) 



Furthermore, the mean number of individuals in the system was found to be 
E [#"(*)] =-E„,i (Af), t>0, 1/6(0,1], (2.3) 

where 

oo r 

f^ r ( ar + P) 

is the generalized Mittag-Leffler function. 

Let W (t) be the number of individuals in a fractional linear birth process 
or fractional Yule or Yule-Furry process (fYp) up to the time t > 0. The state 
probabilities p%(t) = Pr{W(i) = | W(0) = 1} solve the following Cauchy 
problem: 

f &PkV) = - Xk pm + A(fe - iK_i(t), * > i, 

'l, fc = l, (2.4) 
0, fc>l, 
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which is also a fractional generalization of (1.1 1. The fractional derivative 
involved in (2.4) is now the Caputo operator, and is defined as 



'&f(t) = T^tilB&ds, 1/6(0,1), 
/'(<), v =1. 



(2.5) 



Moreover, the Riemann-Liouville (2.2 1 and the Caputo (12. 51) fractional 



derivatives are linked together by the following relation (see Kilbas et al. 
pOOo] , page 91): 



-/w = -/(*)- 



(0,1). 



From (2.6), it is easy to see that both fractional derivatives coincide when 
/(0) = for each k > 1. The solution to the Cauchy problem (2.4) is 

k 



fe- 1 
i - 1 



(-lJ'-'^^-Aif), fc>i, ve{o,i}. 



(2.7) 



Note that the mean number of individuals E[W(t)] in the fractional Yule 
process is the same as (2.3), and the variance can be calculated as 

War («TP(t)) = 2E U>1 (2\t v ) - E uA {\t") - E 2 ^ (At") 

From here on, we emphasize that the fractional derivative operation is per- 
formed in Caputo's sense. 



3 Stretched Yule process with random rates and related 
representations 

In this section, we present some relevant and interesting representations of the 
fractional Yule process (fYp). We start by proving a subordination relation 
that links the fractional Yule process with its classical counterpart. 

Theorem 3.1 Let W(i) be the number of individuals in a fractional Yule 
process at time t > 0. Then the following equality in distribution holds: 

m»(t)lm(T 2u (t)) 7 (3.i) 

where 91(£) is a classical Yule process, v € (0, 1], and Ti v (t) is a random time 
whose distribution coincides with the solution to the following Cauchy problem 

§^g(x, t) = jjgsgix, t), x > 0, 

|:9(M)L = 0, (3-2) 
g(x,0) = S(x), 

with the initial condition g t (x,Q) = 0, when 1/2 < v < 1. 
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Proof Let G u (u,t), t > 0, \u\ < 1, be the probability generating function of the 
fractional Yule process. To prove (3.1), it is sufficient to observe that 

e- zt G v {u,t)dt 

r oo oo 



r _ _ 
/ e-^Ys^pUWt 

J ° k=l 



k - 1 



fe=i i=i 

oo k 



e*<e t-')- 1 ? 



A/ 



e -,(A*+*") ds 



- 1 e - sA '^- 1 e" 



e/eC:;)'-^- 1 / 

/.oo 00 k /h l\ 

X E u *E(*_i)(- 1 ) , " le " A, "/ e--Pr{T 2 ,(t)Gd S }* 

/>oo 

E u k / Pr{*TC(» = fc} Pr{T 2v {t) G ds} dt 

OO 

^« fe Pr{(n(T 2l/ (i))=fc} 



,fc=i 



di. □ 



Remark 3.1 ./Vote i/iat, the solution to (3.2), also solves the fractional differ- 
ential equation 



d v d 
—g( X ,t) = --g( X ,t). 



(3.3) 



Remark 3.2 In the proof of Theorem \3.1\ we used the Laplace transform of 
Pr{T2 V (t) G ds} which is 

/•OO 

/ e' zt Pr{T 2u {i) G ds} = z^e' 8 *" ds, s > 0. 
Jo 

In the next Theorem, we derive a random-rate representation of the frac- 
tional Yule process using the preceding subordination relation. 

Theorem 3.2 (Representation A) Lett > andv G (0, 1]. Then the following 
equality in distribution holds: 



W{t) = m s (f), 



(3.4) 
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where D r ts(i I/ ) is a classical linear birth process with random rate AS evalu- 
ated att v ,Sisa Wright- distributed random variable with probability density 
function W- Vt \- v {—£,) in (1.2). 



Proof To prove equality (3.4 1, we use the subordination relation (3.1) as fol- 
lows: 



Pr{m v (t) = k | 9T(0) = 1} 

poo 

= / Pr{%\{s) = k | 01(0) = 1} Pr{T 2 „{t) e ds} 
Jo 

' k fh — 1 \ 

E ( / _ l ) (- 1 't 1 ^ xlst ^ w ~^(- t ~ Us ) ds 
i=i ^ ' 

1=1 ^ ' 

pOO 

/ Pr{*k<f) = k I 3^(0) = 1} W- v ,i-„(-t)dS, 
Jo 



(3.5) 



and this leads to (3.4 1. □ 



Note that in the second step of formula (3.5 1, we used the explicit form of 



the solution to the fractional diffusion equation (3.2) which is (see Podlubny 
[(19991 ), formula (4.22), page 142) 

PT{T 2v (t) € ds} = r' J W^ lJ . 1 ^ lJ {-t- v s)ds, s > 0. 



Remark 3.3 As noted above, representation (3.4) holds for the one- dimensional 
state probability distribution p%(t), t > 0, k > 1. This, however is sufficient in 
the sense that the process 0Ts(£ w ) has distribution that solves (2.4). 

We now prove a further interesting representation of the fractional Yule 
process in terms of a specific mixed non-homogeneous Poisson process. 
Starting from the second-to-last step of formula (3.5), we obtain 

k 



E 



k-l 
l-l 



^rf-i e -mt" w _ vil _ v{ _ m 



rco k—1 



- 1 



k-l 



Recalling the identity 



c x r dx 



- (r+1) r!, r € N, 91(a) > 0, 
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we get 

pm - 



JO 



1] 



fc-i 



(fc-1) 



oo ,oo „-« e A *« -1 fc-1 f A£t" 



1] 



^0 

oo />oo e 
o Jo 



(fc-1) 



W- V ,x- V (-Z)dud£ (3.6) 
fc-i 



- /p" wA£e A<;s (2s 



fc-i 



(fc-1)! 



Thus, we have obtained a representation in terms of a mixed non-homogeneous 
Poisson process with intensity function 

\{t) = (2\5e xst , t > 0, 

where the distribution of Q is negative-exponential with mean equal to 1, and 
S has probability density function (1.2). Note that the random variable .!?, 
conditional on S — £, is such that 



o. 



as t — v co (see e.g. |Keiding (1974[ ), [Waugh (1970| ), |Harris (2002[ )). 

Remark 3.4 A simple change of variable also allows us to obtain a represen- 
tation in terms of a mixed non-homogeneous Poisson process evaluated at the 
random time Ti v (t\ t > 0. From the second step of formula (3.6), we have 



- W [e*'-l] w fc-l [ e A, _ 



JO 

oo poo 



fc-1 



(fc - 1)! 



-e — 
t 



e- u W- v ,i- v (-€)du d£ 
-W-„,i_„ (-~) dsdu. 



Consider a non-homogeneous Poisson process N(t) with intensity function 
A(t) = f2Ae xt . Then the state probabilities of the fractional Yule process can 
be written as 

/>oo />OG 

P%(t)= e- u Pr{N(s) = k-l}Pr{T 2 „(t)eds}duj (3.7) 
Jo Jo 

= E Q N(T 2l/ (t)), 

In addition, the subordinated non-homogeneous Poisson process N(T2„(i)) 
conditioned on ft = lo could be interesting as the fractional homogeneous 
Poisson process admits a similar representation ( Beghin and Orsingher, 2009[ ). 

Let q%(t) be the state probabilities of N(T2„(t)), i.e., 



g£(i)=Pr{N(T 2l/ (*)) = fc-l}, t>0, fc>l. 



s 
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Then 



oo -w|e* s -l|, ,fc-l \„\s 



fc-1 



(fc-1)! 



-Pr{T 2l/ (t)€ds}. (3.8) 



Applying the Laplace transform to (3.8), we have 



e-*Y k (t)dt -- 
o jo 



oo -w[e AS -l|, ,fc-l f„As 



fc-1 



(fc-1)! 



00 ^e-" 6 'cj^ 1 [l - e Xs ] 
/0 (fc-1)! 

and by taking into account the relations 

00 / , ,\i„AsZ 



fc-1 



z v ~ x e~ sz ds 



(-l) k - 1 ^- 1 e-" r ds, 



E 

5-1 

E 



n 



3=0 



we arrive at the equality 



e- zt ql{t)dt 



(3.9) 



jf ^ (-I)*- 1 -*- 1 E E ( fc ') (-D^v-e-^ 

/» oo 

(-1)V _1 / e- s ^-^ l+ ^ds 



i=o 3=0 

i—n j = Q \ J 



, , fc _ , ~ (- W )'^/fe-l 



■(-^e^eT'tMmf : 



fc-1 



(fc-1) 



1=0 j=0 v 



v-1 



z v -\{l+j)' 



Applying the inverse Laplace transform to equation (3.9), we obtain the ex- 
plicit expression of the state probabilities as 

m = E E ( • Z c-iy- 1 ^ ^ (i + i) n , * > i. 

(3.10) 



Remark 3.5 From equation (3.10), it is straightforward to obtain the classical 



form of the state probabilities of the (conditional) non-homogeneous Poisson 
process (v = 1) with intensity function A(t) = u;Ae Ai , t > 0. 
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We introduce a definition and a lemma below which will be helpful in 
transforming FYp into a non-homogeneous Poisson process with rate 1. In 
order to do so, we present here the standard definition, by means of a Mellin- 
Barnes type integral, of the so-called Fox function: 



Tjm, 71 



where 



(ai,Ai),...,(a p ,A p ) 
(b 1 ,B 1 ),...,(b p ,B p ) 



1 

2?ri 



c-Hoo 



9(z)x~ z dz, x^O, (3.11) 



9(z) = 



1 — aj — AjZ 



} 



(3.12) 



Each empty product is interpreted as unity. For more information on Fox 
functions we refer to Mathai et al. (2010). 

Definition 3.1 Let X^i) be a random time process whose one-dimensional 
distribution is given by 



Pr{1 v {t) e ds} = h(t, s)ds = r^H^l 



(l-i/M/f) 
(o,i) 



ds, 



where t > 0, s > 0, V € (0, 1]. Furthermore, h(t, s) has Mellin transform 

r(rj) 2=i 



s T '- 1 h(t,s)ds = 



(3.13) 



Lemma 3.1 Let y\ u (t) be a fractional Yule process with rate A > and t > 0. 
Then the process 91" (2^ (t)) is a classical Yule process with rate X. 

Proof Define G v (u,i) and G(u,t), t>0, \u\ < 1 as the probability generating 
functions of fYp and the classical Yule process, respectively. Then 



G u (u, s)h(t, s)ds 



r°° 00 k /k — 1 \ 

/ E u " E ■ I (-ly^E^i-xjsnnt, S )d S . 

Jo k=l j=l \ j L s 



Ln the following we use the Mellin-Barnes representation of the Mittag-Leffier 
function 



E v ,i(x) 



1 

2~7ri 



c+ioo 



r(z)rq - z) 
r(i - vz) 



(-x)~ z dz, v > 0, x ^ 



(see Kilbas et al. (2006), page 41, formula (1.8.14))- Note that when v = 1 we 
retrieve the Mellin-Barnes representation of the exponential function 



1 

2tti 



c+ioo 



r{z){-x)- z dz, x^0. 

oo 

(see Paris and Kaminski (200l] l, page 89, formula (3.3.2)). 



(3.14) 
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We obtain 
G v (u, s)h(t, s)ds 



k=l 3 = 1 



k-\\ i-iy- 1 f c+to ° r{z)r{\ - z ) 
^ r(i - vz) 



2iri 



(AjT 



h(t,s) 



ds dz. 



Applying formula (3.13) ; we can write 



G u (u, s)h(t, s)ds 



(3.15) 



= E- fc E 



fc-rw-i)*- 1 r +io ° r(z)r(i-z) .._ z r{i-vz) 



fc=l 3=1 

oo A: 

k=i j=i 

oo ft 

E" fc E 

fe=i j=i 



j — 1 J 2ni 



r(i - z/z) 



(AjT 



r(i-z) 



fc-iv-iy- 1 /• c+io ° 



j — 1 J 2iri 
k-1 



r(z)(\jt)- z dz 



i-i 



J2u k e~ M [l-e~ M ] 



k-l 



k=l 



G(u,t). □ 



Remark 3.6 ./Vote f/ia£ if is straightforward to generalize Lemma 3.1 to the 
more general (non-linear) case. 



Remark 3.7 Letting u — 1 in (3.15), we have 



00 />oo 00 

£ / p%(s)h(t,s)ds = Y,Pk(t) 

k=l J ° k=l 
poo 

^ / h(t, s)ds = 1. 
j 



Theorem 3.3 Consider a fractional Yule process W(t) ttrai/i 6ir£/i raie A > 0, 
i > 0, and ^ € (0, 1]. XTien the random time-changed process 



•r I ^iog 



r? 



ftas one- dimensional distribution which coincides with that of a non-homogeneous 
Poisson process M(t) with rate 1. 



Proof It readily follows from (3.7), Lemma 3.1 and Theorem 1 of Kendall 
JWM). 
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4 Wait and sojourn time distributions 

We now show that the sojourn or inter-birth time of fYp follows the Mittag- 
Leffier distribution. Let T" , i > 1, denote the time between the (i — l)th and 
ith birth. This means that T? is the time it takes for the population size to 
grow from i to i + 1. More specifically, we will show that the sojourn times 
T^'s are independent and Tf is distributed 

f T „(t) = i\t u - 1 E u . v (-i\t u ), i > 1. (4.1) 

Recall that when v = 1, the inter-birth times TVs of the Yp are independent 
and Ti is exponentially distributed with rate i\ i > 1. Moreover, the waiting 
or birth time distribution for the pure linear birth process {v = 1) satisfies the 
following two equalities: 

Pr(2U i = T 1 + ---+T j <t) = Pr((R(t) > j + l|0T(0) = 1) 

and 

p 3 (t) = Pr(2Uj_i < t) - Pr(2Hj < t). 

Let =T? +T% -\ h 7J be the waiting time of the jth birth of the 

fYp. We now show that the preceding two equations hold true as well for the 
fractional or general case (0 < v < 1), i.e., 

Pr^ < t) = Pr(9T(f) > j + l|0T(0) = 1), j>l, (4.2a) 

and 

p"Jt) = Pr(W l ;_ 1 < i) - Pr(2Ir" < t). (4.2b) 



Using (2.7 1, we obtain 



Pr(9T(i) >j + l|DT(0) = 1) = Pr(W(t)=Jfc|W(0) = l) (4.3) 



k=j+i 

3 



1 - Pr(W(t) = k| W(0) = 1) 



fe=i 

3 k 



J- 

This implies that the jth waiting time 2Uj has distribution 

3 k 

I ■ I 



%W = EE(! IlH^W^Mn, t>o, i/€(o,i]. 



fc=l 1=1 

Integrating the preceding equation, we get 

°° 3 k 

j k=i i=i v 7 
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j fc-1 

EE 

k=l 1=0 



k - 1 



(-!)' = £(1 -l) k - l = l. 



k=l 



The non- negativity of fw<(t) follows from the non- negativity of p£(i) (see 



Orsingher and Polito (2010 )), and the last line of (4.3 ) is a monotone increasing 



function of t. To see this, we can write 



j k 



i-ee(J_i 1 )c- i ) , - i ^(-^) 



k=l 1 = 1 




i-E^w 



k=l 
j k 



i-EE 



fe=i i=i 

j k 



EE 

fe=i i=i 



k-l 
I - 1 

fc-1 



(-lj'-^rCT, > t) 



Z - 1 



(-l^Pr^ < t). 



Indeed, f<xo v (t) is a probability density function. Note also that /gjjv (t) has the 
following integral representation: 



OO 3 k r, , 

e" f EE 

fe=l !=1 



Z - 1 



(-lJ'-^CiAt/Odf, 



where 3(77) = sin(^7r)/[7r(ry l/ + ry _!y + 2 cos(ptt))] (see Repin and Saichev (2000)). 
We now show that if the sojourn times are distributed as in (4.1 1, the cumula- 
tive distribution function Pr(2UJ < t) of the waiting or birth time equals the 
right-hand side of (4.2a). When j = 1, we get 

Pr(2H^ < t) = Pr(T^ < t) = 1 - E v>1 (-Xt u ) = 1 - 

In the succeeding calculations, we use the following identities (see page 26 of 



Podlubny (1999)): 



E Uj i(-jX(t-u) u )u"- 1 E 1/tl/ (-Xlu ,/ )du = 



jE u>v+ i{-j\t») - IE V , V+1 (-I\r)„ 



3-1 



and 



Now, 



l<j. 



Pr(22J£ < t) = f Pr{T? + T% < t\T? = u}dF Tr 
Jo 



(«) 



[1 - E vA (-2X(t - uY)}Xu v - x E v ^{-Xu v )du 
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= 1 - E Vtl {-\t v ) - [2Af E v>v+1 (-2Xt u ) - t u E v>1/+1 (-Xt v )] 
= 1 - E Uil (-Xt u ) - [E u<1 (-Xt u ) - E v>1 (-2Xt v )\ 
= 1 - 2E V>1 (-Xt u ) + E Vtl {-2Xt v ) 

2 

k=l 



and in general, we can show by induction that 
Pr(2D£ < t) 

= f Pr{2D£ < t\W^_ 1 = u}dF w , » 
Jo 

[l-E„ A (-jX(t-uy)]fm, U (u)du 

3-1 k 



[ [1 - E vA {-jX{t - «)")] EE (L !) (-lj'-^AOt*"- 1 ^,^ 

J ° fc=l i=l ^ ' 



-Xlu v )du 



j-i k 



k=l 1=1 
3-1 fc 



EEltJlHl'-'Ii-i.H"')] 



eeC:;)*- 1 '" 1 "/^ 

fe=i /=i v 7 ■ / ° 



(-jA(i - u) v )u v - l E^ u (-Xlu v )du 



3-1 fc 



fe=l i=l 

3-1 fc 



EEt^^-^ti-^i^m 



fc=i /=i 

3-1 k 



l-l I 



J -I 



[E^i-Xlf) - E uA (-Xjt v ) 



k=l 1 = 1 
3-1 k 



k-l\, / j 



. ; £„,i(-AJf)- -— E Ujl (-Xjt v ) 
J -I J -I 



k=l 1=1 
3-1 fe 



fc - 1 



fc=l 2=1 v 7 



X 2-l J 

i 



E Vfl {-Xlt v ) 



' 1 ' E uA (-Xjt"). 



Using the formulas on page 3 of Gradshteyn and Ryzhik (1980), we have 
3-1 k 



k=l 1 = 1 



EE(t> 



I 



3-1 



" ft=2 v 



;=i 
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a -!)!(.? -0! 



because 



Hence, 

Pr(2D£ < t) (4.4) 



1 - EE ft _ 1 1 )(-i)'- 1 7 f 7 ^,i(-A/n - (-lr^.aC-A^) 
fe=i z=i ^ 7 ^ 

i-EEC: 1 1 ) ( - 1) "^ ( - Ain 

*;=i i=i v 7 J 

- (e ({ Z J) (-ij'-^.iC-a/o - E I J) (-i)'-Xi(-A/n 



u=i x 7 2=1 



*;=i i=i v ' x ' 
- (E ({ Z J) (-i) , - 1 ^,i(-Aio - E ( j 7 ') (-i)'- 1 j^^ic-Ain 



u=i x 7 2=1 

i-ei:0: 1 1 )(- i )'- i A e -<-^» 



fc=l 2=1 
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3-1 k 



k=l 1=1 
/ 3 



fc - 1 

Z - 1 



\z=l 

3-1 fc 



j - 1 

3-1 fc 



E^Zjjc-D'-Xi^n-EE 

fc - 1 



l=k 1=1 
3 



=1 -EE ; H)'%(-^)-E 
fe=i i=i ^ ' i=i 

3 

= l-$>fc(*)» l<fc<j ! 



jfe - 1 
Z- 1 

.7 - 1 
Z - 1 



(-IJ'-'j^^xC-AZO 



(-^'-^(-Alt") 



fc=i 



as the second summation (in the preceding equal sign) simply corresponds 
to k = j. Hence, equality (4.2a| is attained. Again, the transition from the 
third equality to the fourth equality above uses formula (0.15.1) on page 3 of 



Gradshteyn and Ryzhik (19801, i.e., 

3-1-1 



3-1 

E 

k=l 



jfe-1 
I - 1 



E 

fe=0 



k + l-1 
l-l 



Notice that when v = 1, we get Pr(2Uj < t) = (1 — e which corresponds 
to the birth time distribution of the classical Yule process. Moreover, equality 



(4.2b) can be straightforwardly evaluated as 

/ 3-1 fc 



Pr(22JJ_ x < t) - Pr(2»J < *) = 1 - EE 

V fc=i (=i 

-fi-EE 



fc=i ;=i 



E 

£=1 



.7-1 
Z - 1 



fc- 1 
I - 1 

fc- 1 
l-l 

\l-l 



{-l) l - l E v , x {-\lt v )j 
{-l) 1 - 1 E V:l {-\lt v ) 



(-iy- i E Utl (-\it v ) 



In addition, the Laplace transform of the probability density fxv {t) is 

iX 



POO 

/ e~ zt 
Jo 



hAt)dt = 



i\ + z v 



This suggests that the distribution (eqn (4.1)) leads to the following known 
mixture or structural representation (see Cahoy et al. (2010 )) of the inter-birth 
times as 

If = V? /v S v , 

where Vi has the exponential distribution with parameter iX, i.e., 

f v .( v ) =iXe- lXv , v>0, (4.5) 



16 



Dexter O. Cahoy, Federico Polito 



and is independent of the positive Levy or (/-stable distributed random vari- 
able S v having the Laplace transform of the density function e~ z . This also 
suggests that the n-th fractional moment of the ith inter-birth time is given 

by 

[lJ (i\)*r(K/v) Bin(7r«/i/)r(l - /s) ' ' 

which further implies that the n-th fractional moment of the jth wait or birth 
time is 

F rorr^l « = 7tT(1 + /s) ^ A /fc - l\ ^ / 1_ 

L J ' J A K r( K /z/)sin(7TK/^)r(i- K ) V- 1 / V" 

where < k < v. 



5 Sample paths of fYp 

From Sections 3 and 4, it is now straightforward to simulate a trajectory of a 
fYp. However, we only propose the two simplest algorithms on how to generate 
a sample path of the fYp as the others follow. In particular, the random-rate 



representation (Representation A, Theorem 3.2 ) yields the algorithm below. 
ALGORITHM 1: 

i) Generate S from the Wright distribution W- u ,i-v(—£), and obtain £. 

ii) Simulate a classical Yule process with birth rate A£. 

iii) Stretch the time scale to t". 

A simpler way to generate a realization of fYp with n births is to directly 
exploit the known birth and/or sojourn time distributions as follows: Generate 



Vi from the exponential distribution in (4.5) with parameter iX, and S u from 
the strictly positive stable distribution with parameter v. 

ALGORITHM 2: 

i) Let i = 1 and 0F(0) = 1. 

ii) Simulate T" = V^'" S v , and let = + T% + • • • + . 

iii) W(2Uf) = % + 1, and % = i + 1. 

iv) Repeat ii-iii for i = 2, . . . , n — 1. 

We now use the algorithms above to highlight some unique properties of 



the fractional Yule process that are related to its true mean given in (2.3). 
Figure |5.1| below shows both Yp and fYp as jump processes of size 1 in the 
time interval (0, 5) with v — 0.5, and A = 1. Using the same set of parameters, 



Figure 5.2 displays sample trajectories of a different/independent fYp and Yp 
which model a binary-split growth process. An important attribute that can 
be directly observed from these two graphs is that on the average, fYp grows 
more rapidly than the classical Yp shortly after it starts. 
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Fig. 5.1 Sample trajectories of the standard Yule process (top) and the fractional Yule 
process (bottom) in the interval (0, 5) with parameters (u, A) = (0.5, 1). 




In addition, a more specific characteristic of fYp is illustrated in Figure 5.3 



The particular realization of fYp below used the parameter values v — 0.25, 
A = 1, and is observed in the time interval (0, 5) . It clearly suggests that fYp 
is more explosive than Yp when v — \ 0. In general, the plots strongly validate 
the plausibility of fYp to model exploding and strictly growing processes. Note 
also that Representation A implies that the interaction between the random 
rate and time stretching of the classical Yule process can rapidly speed up or 
slow down fYp at any given time instance. 
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6 Method-of-Moments (MoM) estimation 

We now propose a method-of-moments estimation procedure for the param- 
eters v and A to make fYp usable in practice. In this procedure, we assume 
that a particular realization or complete history of the process is observed un- 
til the population is n, i.e., there are n births. We then attempt to use all the 
available data from the observed sample path of the fractional Yule process. 

In particular, we use all the available inter-birth or sojourn times of the 
observed sample trajectory of the fractional Yule process. A direct way of es- 
timating the parameters is to use the fractional moment estimators as follows: 
Choose constants n m < v : m = 1, 2, and solve for the estimates A and v using 
the equations 



7lT(l 



n 



\ K ™r(n m /v) six\(iTK m /v)r(l - K m ) 



1,2. 



Another approach is to use the first two integer-order moments of the log- 



transformed sojourn times (see Cahoy et al. (2010)) which are 



ElnPTH^-7, 



and 



Eln [T»] 2 = tt 2 



Zv" 1 



ln(iA) 
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This further suggests that the parameter estimates can be computed using the 
two equations: 

EIUinp7]_-£r=iM*A) 



n vn 

and 



-7, 



n = " {M-6j + n^{— +1 

where 7 = 0.577215664901532 is the Euler-Mascheroni constant. A major 
advantage of this procedure over other moment estimators is that it does not 
require selection of constants a priori to calculate the parameter estimates. 
Note also that the maximum likelihood estimators are more challenging to 
compute due to the required evaluation of the Mittag-Leffler function. 

In addition, we tested our parameter estimation procedure. In doing so, 
we generated 10 random samples of inter-birth times of size 10000 each for 
v = 0.1 + 0.1m, m = 0, . . . , 9 and A = 0.2, 10. For each simulated data set, 
we computed the estimates using the first n observations in the set with n = 
100, 1000, and 10000. The tables below show the simulation results for a single 
run, which further indicate that the proposed procedure performs relatively 
well as the sample sizes increase. Please note that in many applications (e.g., 
internet traffic), the typical number of observations is at least of the order 
of millions. These estimates could also serve as good starting values of an 
iterative estimation procedure. 



Table 6.1 Parameter estimates (v, A) for fYp with v = 0.1(0.1)1 and A = 0.2. 





n = 100 


n = 1000 


n = 10000 


[y 


= 0.1, A 


= 0.2) 


(0.095, 0.198) 


(0.096, 0.185) 


(0.100, 0.205) 


[y 


= 0.2, A 


= 0.2) 


(0.228, 0.249) 


(0.193, 0.189) 


(0.199, 0.193) 


[y 


= 0.3, A 


= 0.2) 


(0.283, 0.185) 


(0.292, 0.193) 


(0.303, 0.228) 


(y 


= 0.4, A 


= 0.2) 


(0.381, 0.178) 


(0.407, 0.218) 


(0.402, 0.209) 


[y 


= 0.5, A 


= 0.2) 


(0.481, 0.212) 


(0.501, 0.197) 


(0.500, 0.197) 


[y 


= 0.6, A 


= 0.2) 


(0.599, 0.211) 


(0.602, 0.186) 


(0.595, 0.186) 


[y 


= 0.7, A 


= 0.2) 


(0.759, 0.257) 


(0.728, 0.250) 


(0.700, 0.198) 


[y 


= 0.8, A 


= 0.2) 


(0.818, 0.220) 


(0.819, 0.229) 


(0.803, 0.204) 


[y 


= 0.9, A 


= 0.2) 


(0.850, 0.193) 


(0.899, 0.211) 


(0.907, 0.215) 


(y 


= 1.0, A 


= 0.2) 


(0.977, 0.183) 


(0.991, 0.199) 


(0.999, 0.202) 



7 Concluding remarks 

We have derived one-dimensional representations of the fractional Yule pro- 
cess, which led to algorithms for simulating its sample paths. These repre- 
sentations are also necessary in understanding the properties of fYp further. 
We have derived the birth and inter-birth or sojourn time distributions, which 
are of Mittag-Leffler type. The structural representation of the random so- 
journ time also led to an algorithm for simulating sample trajectories of the 
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Table 6.2 Parameter estimates (P, A) for fYp with v = 0.1(0.1)1 and A = 10. 





n = 100 


n = 1000 


n = 10000 


[y 


= 0.1, A = 


10) 


(0.107, 13.067) 


(0.101, 10.599) 


(0.101, 10.730) 


[y 


= 0.2, A = 


10) 


(0.203, 10.737) 


(0.206, 12.384) 


(0.201, 10.555) 


[y 


= 0.3, A = 


10) 


(0.299, 11.027) 


(0.297, 9.359) 


(0.295, 8.593) 


{y 


= 0.4, A = 


10) 


(0.391, 7.598) 


(0.396, 8.899) 


(0.397, 9.086) 


[y 


= 0.5, A = 


10) 


(0.517, 10.939) 


(0.509, 11.428) 


(0.501, 10.269) 


[y 


= 0.6, A = 


10) 


(0.630, 11.379) 


(0.586, 8.308) 


(0.597, 9.162) 


[y 


= 0.7, A = 


10) 


(0.716, 12.413) 


(0.699, 10.634) 


(0.710, 11.679) 


[y 


= 0.8, A = 


10) 


(0.782, 8.713) 


(0.786, 8.186) 


(0.804, 10.498) 


[y 


= 0.9, A = 


10) 


(0.919, 11.429) 


(0.899, 9.043) 


(0.897, 9.684) 


[y 


= 1.0, A = 


10) 


(0.969, 8.712) 


(1.000, 10.427) 


(1.001, 10.434) 



fYp. Wc have proposed an estimation procedure using the moments of the 
log-transformed inter-birth times, which performed satisfactorily especially for 
larger sample sizes. 

Although some properties of fYp have already been studied, there are still 
a lot of open problems that need to be figured out. For instance, understanding 
fYp in more depth and the construction of more efficient estimators like the 
maximum likelihood would be worth pursuing in the future. Also, the appli- 
cation of fYp in practice particularly in biology and/or network traffic is still 
in progress. 
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